function [Lambda_hat_AJX, Vn_AJX, beta_J_hat,beta_C_hat,Lambda_hat_AJX_truebeta] ...
    = AJX_HF(dPt, dFt, L_k, B_kl, Delta_n, un, qn, vn, JumpVectorF, JumpVectorX, jump, beta_C_tr, beta_J_tr)

% This code allows nan for missing data, but require qn to be divisor of n
% for speed and better finite sample performance.

% L_k: the input should be a list, with number of jump factors for each
% factor; sum of L_k = H
% B_kl: the input should be a H*2 matrix. Specifying the range of the H
% jump factors.

    K         =      size(dFt,1);
    H         =      sum(L_k);
    M         =      size(dPt,1);
    n         =      size(dPt,2);
    Tn        =      round(n*Delta_n);
    nobs      =      (n/(Tn*252));
    dPt_zero  =      dPt;
    dPt_zero(isnan(dPt)) = 0;

    lambdat_hat       = zeros(K+H,(floor(n/qn)-1));
    lambdat_hat_truebeta  = zeros(K+H,(floor(n/qn)-1));
        Vnt_hat       = zeros(K+H,K+H,(floor(n/qn)-1));
        beta_C_hat    = zeros(M,K,(floor(n/qn)-1));
       
    % Jump regression
    % Use this matrix to filter factor jump to make sure we are considering
    % time period corresponding to each individual stock's live time. 
    dFt_filter_P = zeros(H, n, M);
    for i = 1:M
        dFt_filter_P(:,:,i) = repmat(~isnan(dPt(i,:)), H, 1);
    end
    
    % Step 1 Stack all the factor jump together into dFt_bar(H*n)
    if jump == false
        dFt_cont_ix = fn_cont_ix(un, 0.47, dFt, Delta_n, Tn*252, nobs);
    else
        dFt_cont_ix = [JumpVectorF == 0];
    end
    dFt_jp_ix = ones(K,n) - dFt_cont_ix;
    
    dFt_bar = zeros(H,n);
    id = zeros(H,n);
    
        
    for i = 1:K
        for j = 1:L_k(i)
            m = sum(L_k(1:i-1));
            id(m+j, :) = dFt_jp_ix(i,:) & dFt(i,:)>min(B_kl(m+j,:)) & dFt(i,:)<max(B_kl(m+j,:)); 
        end
    end
    
    for i = 1:K
        for j = 1:L_k(i)
            m = sum(L_k(1:i-1));
            dFt_bar(m+j,:) = dFt(i,:);
        end
    end
    
    dFt_bar = dFt_bar.*id;
    beta_J_hat = zeros(M,H);
    for i = 1:M
        if min(eig((dFt_bar.*dFt_filter_P(:,:,i))*(dFt_bar.*dFt_filter_P(:,:,i))')) > 1/vn
            beta_J_hat(i, :) = dPt_zero(i,:)*(dFt_bar') *...
                            inv((dFt_bar.*dFt_filter_P(:,:,i))*(dFt_bar.*dFt_filter_P(:,:,i))');
        else
            beta_J_hat(i, :) = 0;
        end
    end
        
    if jump == true
        dFt_true = dFt.* ~JumpVectorF;
        dPt_zero_true = dPt_zero .* ~JumpVectorX;
    end
    % beta_J_hat = beta_J;
    for j = 0:(floor(n/qn)-2)
        
        
        Ij = (j*qn+1):((j+1)*qn);
        Ijp = ((j+1)*qn+1):((j+2)*qn);
        if jump == true
            dFt_j = dFt_true(:,Ij);
            dPt_j = dPt_zero_true(:,Ij);
        else
            dFt_j = dFt(:, Ij);
            dPt_j = dPt_zero(:,Ij);
        end
        
        % Deal with NaN
        for i = 1:M
            if sum(isnan(dPt(i,Ij))) > 0
                beta_J_hat_j(i,:) = zeros(1,H);
            else
                beta_J_hat_j(i,:) = beta_J_hat(i,:);
            end
        end
        
        % Filter the jump 
        if jump == false
            dFt_cont_ix_j = fn_cont_ix(un, 0.47, dFt_j, Delta_n, qn/nobs, nobs);
            dPt_cont_ix_j = fn_cont_ix(un, 0.47, dPt_j, Delta_n, qn/nobs, nobs);
        else
            dFt_cont_ix_j = ones(size(dFt_j));
            dPt_cont_ix_j = ones(size(dPt_j));
        end
        
        % Retain only continuous movements in factors
        dFt_j = dFt_j.*dFt_cont_ix_j;
        dPt_j = dPt_j.*dPt_cont_ix_j;
        
        cF_t    = 1/(qn*Delta_n) * dFt_j*(dFt_j');
        gamma_t = 1/(qn*Delta_n) * (dPt_j*dFt_j');
        
        if min(eig(cF_t)) > 1/vn
            beta_C_hat_j = gamma_t * inv(cF_t);
        else
            beta_C_hat_j = zeros(M,K);
        end
        
        for i = 1:M
            if sum(isnan(dPt(i,Ij))) > 0
                beta_C_hat_j(i,:) = zeros(1,K);
            end
        end
        
        beta_hat = [beta_C_hat_j, beta_J_hat_j];
        b_hat = beta_hat' * beta_hat;
        
    %    if min(eig(b_hat)) > 1/vn
            eta_hat = pinv(b_hat) * beta_hat'; % If there is no jump beta, sets its risk premia to 0, but other part can still be estimated.
         %  else
         %   eta_hat = zeros(K+H,M);
        %end
        lambdat_hat(:,j+1) = eta_hat * sum(dPt_zero(:,Ijp),2); 
        Vnt_hat(:,:,j+1) = eta_hat * sum(dPt_zero(:,Ijp),2)*(sum(dPt_zero(:,Ijp),2)') * (eta_hat');

        beta_C_hat(:,:,j+1) =  beta_C_hat_j;
        
        %%%%%
        beta_C_tr_j = beta_C_tr(:,:,j*qn+1);
        beta_tr_j = [beta_C_tr_j, beta_J_tr];
        
        lambdat_hat_truebeta(:,j+1) =  inv(beta_tr_j'*beta_tr_j) * beta_tr_j' * sum(dPt_zero(:,Ijp),2);
        
        
        
    end
    
    Lambda_hat_AJX_truebeta = sum(lambdat_hat_truebeta,2)/(floor(n/qn)*qn*Delta_n); 
    Lambda_hat_AJX = sum(lambdat_hat,2)/(floor(n/qn)*qn*Delta_n); 
    Vn_AJX         = sum(Vnt_hat,3)/(floor(n/qn)*qn*Delta_n)^2;
end